import requests
from datetime import datetime, timedelta, timezone
import os, json
import matplotlib.pyplot as plt
import matplotlib.dates as mdates
import numpy as np

from astropy.coordinates import SkyCoord, AltAz, EarthLocation, ICRS
from astropy.utils.iers import IERS_Auto, IERS_A
from astropy.time import Time
from astropy import units as u
# from astropy_healpix import HEALPix

'''
finals2000A.all下载失败：可以使用如下3种方法
'''
#1， 设置 IERS_Auto，使其尽量使用本地缓存
# iers = IERS_Auto.open()

#2， 读取本地旧数据
iers_old = IERS_A.read("public/svomSummary/python/finals2000A_20250321.all")
from astropy.utils.iers import conf # 设置 IERS 表，使 astropy 使用这个数据
conf.iers_auto_url = None  # 禁止自动下载新数据
conf.iers_table = iers_old

#3，禁用 IERS 自动更新（使用默认值）
# from astropy.utils.iers import conf
# conf.auto_max_age = None  # 禁用 IERS 自动更新

def format_date(date):
    return date.strftime('%Y-%m-%dT%H:%M:%S')

def get_current_date():
    return datetime.now()

def get_date_before_days(days, dateStr=""):
    if len(dateStr) > 0:
        dateObj = datetime.strptime(dateStr, "%Y-%m-%d")
    else:
        dateObj = get_current_date()
    
    return dateObj - timedelta(days=days)

def get_date_after_days(days, dateStr=""):
    if len(dateStr) > 0:
        dateObj = datetime.strptime(dateStr, "%Y-%m-%d")
    else:
        dateObj = get_current_date()
    return dateObj + timedelta(days=days)

def request_url(root_url, dateStr, begin_time_param_name="begTime"):
    two_days_before = get_date_before_days(2, dateStr)
    one_day_after = get_date_after_days(1, dateStr)

    two_days_before_str = format_date(two_days_before)
    one_day_after_str = format_date(one_day_after)

    url = f"{root_url}?{begin_time_param_name}={two_days_before_str}&endTime={one_day_after_str}"
    # print(url)

    data_list = []
    try:
        response = requests.get(url, timeout=10)
        response.raise_for_status()
        data = response.json()
        
        if data.get('errorCode') == 0 or data.get('errCode') == 0:
            data_list = data.get('data', [])
    except requests.RequestException as e:
        print(f"请求失败: {e}")
    
    return data_list


def getSvomWorkplan(saveDir, dateStr):
    
    url = "https://www.svom.cn/svomGP/obs-corr-info/getFmtdDurObsPlans"
    # svom_url = "https://www.svom.cn/svomGP/obs-corr-info/getFmtdDurObsPlans?begTime=2025-02-26T00:00:00&endTime=2025-02-27T00:00:00"
    # svom_url = "https://www.svom.cn/svomGP/obs-corr-info//getDurObsPlans?begTime=2025-02-26T00:00:00&endTime=2025-02-27T00:00:00"

    svom_data_file = f"{saveDir}/svom_data.json"
            
    if not os.path.exists(svom_data_file):
        print('fetching data from SVOM...')
        data_list = request_url(url, dateStr, 'beginTime')
        with open(svom_data_file, "w") as f:
            json.dump(data_list, f, indent=4)
    else:
        with open(svom_data_file, "r") as f:
            data_list = json.load(f)
    
    priority_order = {'TOO-EX': 3, 'TOO-NOM': 2, 'OTHER': 1, 'GP': 0}
    known_categories = ['GP', 'TOO-EX', 'TOO-NOM']
    
    tdata = []
    ttcategorys = []
    for element in data_list:
        tcategory = element.get('cls', '')
        start_date_str = element.get('beginTime', '') + "Z"
        end_date_str = element.get('endTime', '') + "Z"
        if not start_date_str or not end_date_str or len(start_date_str) == 0 or len(end_date_str) == 0:
            continue
        start_date = datetime.fromisoformat(start_date_str[:-1])
        end_date = datetime.fromisoformat(end_date_str[:-1])
        
        info = element.get("info", {})
        gsp_id = info.get("gspId", None)
        ra = info.get("rightAscension", None)
        dec = info.get("declination", None)

        tcategory = tcategory.replace("-work-plan", "").upper()
        if tcategory not in known_categories:
            tcategory = 'OTHER'
        if tcategory not in ttcategorys:
            ttcategorys.append(tcategory)
        tdata.append({
            'priority': priority_order.get(tcategory, -1),
            'category': tcategory,
            'eventDate': start_date,
            "start": start_date,
            "end": end_date,
            "ra": ra,
            "dec": dec
        })
        
    # print(ttcategorys)
    
    return tdata
def filter_overlapping_events(events):
    # 先按照优先级从高到低排序，确保高优先级事件先加入
    events.sort(key=lambda x: (-x['priority'], x['start']))
    filtered_data = []

    for item in events:
        new_segments = [(item['start'], item['end'])]

        for existing in filtered_data:
            if existing['end'] <= item['start'] or existing['start'] >= item['end']:
                # 没有重叠，跳过
                continue

            if existing['priority'] >= item['priority']:
                # 当前 item 是低优先级，需要去掉重叠部分
                updated_segments = []
                for start, end in new_segments:
                    if start < existing['start']:
                        updated_segments.append((start, existing['start']))
                    if end > existing['end']:
                        updated_segments.append((existing['end'], end))
                new_segments = updated_segments

        # 添加处理后的非重叠部分
        for start, end in new_segments:
            if start < end:  # 避免无效时间段
                new_item = item.copy()
                new_item['start'] = start
                new_item['end'] = end
                filtered_data.append(new_item)

    return filtered_data

def find_overlapping_events(events, threshold=60):
    overlapping_pairs = []
    events.sort(key=lambda x: x['start'])
    
    for i in range(len(events)):
        for j in range(i + 1, len(events)):
            if events[j]['start'] >= events[i]['end']:
                break  # 后续事件不会重叠，跳出内循环
            overlap_start = max(events[i]['start'], events[j]['start'])
            overlap_end = min(events[i]['end'], events[j]['end'])
            
            if (overlap_end - overlap_start).total_seconds() > threshold:
                overlapping_pairs.append((
                    (events[i]['category'], events[i]['start'], events[i]['end']),
                    (events[j]['category'], events[j]['start'], events[j]['end'])
                ))
    
    return overlapping_pairs

def adjust_event_intervals(events):
    events.sort(key=lambda x: (x['start'], -x['priority']))
    min_interval = timedelta(minutes=1)
    events = [event.copy() for event in events]  # 复制数据，防止修改原列表
    
    for i in range(len(events) - 1):
        current_event = events[i]
        next_event = events[i + 1]
        
        if next_event['start'] - current_event['end'] < min_interval:
            if next_event['priority'] < current_event['priority']:
                next_event['start'] = current_event['end'] + min_interval
            else:
                current_event['end'] = next_event['start'] - min_interval
    
    return events

def getProcessedSvomWorkplan(saveDir, dateStr):
    workplan_data = getSvomWorkplan(saveDir, dateStr)
    filtered_data = filter_overlapping_events(workplan_data)
    adjusted_data = adjust_event_intervals(filtered_data)
    
    overlapping_pairs = find_overlapping_events(adjusted_data, threshold=60)
    for pair in overlapping_pairs:
        print(f"Warning: Overlapping Pair: {pair[0][0]} and {pair[1][0]}, Start Time: {pair[0][1]}, End Time: {pair[0][2]}, Start Time: {pair[1][1]}, End Time: {pair[1][2]}")
        
    return adjusted_data

    
def spherical_cap_area2(radius):
    """ 计算球面区域面积 """
    return 2 * np.pi * (1 - np.cos(radius.to(u.rad).value))

def spherical_cap_area(radius):
    """计算球冠面积"""
    return 2 * np.pi * (1 - np.cos(radius.to(u.rad).value)) * u.steradian

'''
基于HEALPix采样模拟地平30度以上的区域，计算重叠面积
'''
def compute_overlap_xinglong_healpix(svom_obs):
    """计算地平30度以上与SVOM指向的重叠比例（相对 SVOM 总视场）"""
    
    obs_location = EarthLocation.from_geodetic(117.575, 40.393, 960)  # 兴隆观测站
    
    overlap_ratios = []
    svom_field_radius = 90 * u.deg  # SVOM 视场半径为 90 度
    min_altitude = 30 * u.deg  # 兴隆可见区域下限
    cap_above_horizon_radius = 60 * u.deg
    
    svom_field_area = spherical_cap_area(svom_field_radius)  # SVOM 视场总面积
    visible_area = spherical_cap_area(cap_above_horizon_radius)  # 地平 30 度以上区域的面积
    
    # 时间范围调整为 UTC 9:00 - 23:00
    utcTime = datetime.now(timezone.utc) - timedelta(days=1)
    start_time = utcTime.replace(hour=10, minute=30, second=0, microsecond=0)
    end_time = utcTime.replace(hour=22, minute=30, second=0, microsecond=0)
    
    selSvomObs = []
    for tobs in svom_obs:
        start, end, center, _ = tobs
        if start.tzinfo is None:
            start = start.replace(tzinfo=timezone.utc)
        if end.tzinfo is None:
            end = end.replace(tzinfo=timezone.utc)
        
        #起止时间是否包含在SVOM的一次观测时间内
        if start <= start_time <= end or start <= end_time <= end: 
            selSvomObs.append(tobs)
        
        #SVOM的一次观测时间是否包含在时间范围内
        if start_time <= start <= end_time and start_time <= end <= end_time:
            selSvomObs.append(tobs)
    
    # for start, end, center, _ in selSvomObs:
    #     print(start, end, center.ra.deg, center.dec.deg)
    
    # 每 10 分钟采样一次
    tsteps = int((end_time - start_time).total_seconds() / 600)
    times = [start_time + timedelta(minutes=i * 10) for i in range(tsteps)]

    # 使用 HEALPix 网格提高采样精度
    nside = 64  # 可根据精度调整
    hp = HEALPix(nside=nside)
    ra, dec = hp.healpix_to_lonlat(np.arange(hp.npix))
    coords = SkyCoord(ra, dec, unit='deg')

    for time in times:
        astro_time = Time(time)
        frame = AltAz(obstime=astro_time, location=obs_location)
        altaz = coords.transform_to(frame)

        visible = altaz.alt > min_altitude  # 筛选出地平 30 度以上的区域
        
        # 检查当前时间是否有 SVOM 观测
        for start, end, center, _ in svom_obs:
            # 将 start 和 end 转换为有时区信息的 datetime 对象
            if start.tzinfo is None:
                start = start.replace(tzinfo=timezone.utc)
            if end.tzinfo is None:
                end = end.replace(tzinfo=timezone.utc)
            
            if start <= time <= end:
                # 计算 SVOM 视场与兴隆可见区域的重叠
                sep = coords[visible].separation(center)
                
                overlap = np.sum(sep < svom_field_radius) / len(visible)
                overlap_area = overlap * visible_area
                
                # 计算相对重叠比例
                overlap_ratio = overlap_area / visible_area
                overlap_ratios.append((time.isoformat(), overlap_ratio))
                break  # 同一时间段只有一个条目，直接跳出循环
    
    return overlap_ratios


def compute_svom_overlay_local_zenith(svom_center, curTime, location):
    
    zenith_radius = 60 * u.deg
    svom_radius = 90 * u.deg
    
    # 转换为astropy Time对象
    curTime = Time(curTime)
    
    # 计算当前天顶的赤道坐标
    altaz = AltAz(alt=90*u.deg, az=0*u.deg, location=location, obstime=curTime)
    zenith_coord = altaz.transform_to(ICRS())
    
    # 计算天顶圆与SVOM圆的角距离
    separation = zenith_coord.separation(svom_center)
    d_deg = separation.to(u.deg).value
    r1_deg = zenith_radius.to(u.deg).value
    r2_deg = svom_radius.to(u.deg).value
    
    # 计算重叠面积
    if d_deg >= r1_deg + r2_deg:
        overlap_area = 0.0
    elif d_deg <= abs(r1_deg - r2_deg):
        overlap_area = np.pi * (min(r1_deg, r2_deg) ** 2)
    else:
        # 转换为弧度计算
        d_rad = np.deg2rad(d_deg)
        r1_rad = np.deg2rad(r1_deg)
        r2_rad = np.deg2rad(r2_deg)
        
        part1 = r1_rad**2 * np.arccos((d_rad**2 + r1_rad**2 - r2_rad**2) / (2 * d_rad * r1_rad))
        part2 = r2_rad**2 * np.arccos((d_rad**2 + r2_rad**2 - r1_rad**2) / (2 * d_rad * r2_rad))
        part3 = 0.5 * np.sqrt((-d_rad + r1_rad + r2_rad) * (d_rad + r1_rad - r2_rad) * 
                                (d_rad - r1_rad + r2_rad) * (d_rad + r1_rad + r2_rad))
        overlap_area_rad_sq = part1 + part2 - part3
        overlap_area = overlap_area_rad_sq * (180.0 / np.pi) ** 2  # 转换为平方度
    
    # 计算比例
    total_zenith_area = np.pi * (r1_deg ** 2)
    ratio = overlap_area / total_zenith_area if total_zenith_area > 0 else 0.0
    
    return ratio
    
'''
用天顶为中心半径60度的圆代表地平30度以上的天区域，计算重叠面积
'''
def compute_overlap_xinglong(svom_obs, interval_minutes=1):
    """
    计算天顶圆与SVOM观测区域的重叠比例随时间的变化。
    
    参数:
    svom_obs (list): SVOM观测列表，每个元素为(start, end, center, category)，
                    其中start/end为datetime，center为SkyCoord。
    location (EarthLocation): 观测地点的地理位置。
    start_date (datetime): 观测的日期，时间部分将被忽略，从当天09:00开始。
    interval_minutes (int): 采样间隔分钟数，默认为10分钟。
    
    返回:
    list: 每个时间点的重叠比例，包含时间及比例。
    """
    location = EarthLocation.from_geodetic(117.575, 40.393, 960)

    interval = timedelta(minutes=interval_minutes)
    
    # 生成时间序列：从当天09:00到23:00
    # utcTime = datetime.now(timezone.utc) - timedelta(days=1)
    utcTime = datetime.now() - timedelta(days=1)
    start_time = utcTime.replace(hour=10, minute=30, second=0, microsecond=0)
    end_time = utcTime.replace(hour=22, minute=30, second=0, microsecond=0)
    
    current_time = start_time
    times = []
    while current_time <= end_time:
        times.append(current_time)
        current_time += interval
    
    results = []
    
    for curTime in times:
        
        # 查找当前时间的SVOM观测
        svom_entry = None
        for entry in svom_obs:
            obs_start, obs_end, svom_center, category = entry
            if obs_start <= curTime <= obs_end:
                svom_entry = entry
                break
        if svom_entry is None:
            # results.append((time.isoformat(),0.0))
            continue
        
        ratio = compute_svom_overlay_local_zenith(svom_center, curTime, location)
        results.append((curTime.isoformat(), ratio))
    
    return results


def plot_workplan2():
    
    # utcTimeStr = tools.getCurDateStr()
    utcTimeStr = "2025-02-26"
    dataDir = f"public/svomSummary/data/{utcTimeStr}"
    
    workplan_data = getSvomWorkplan(dataDir)
    filtered_data = filter_overlapping_events(workplan_data)
    # filtered_data = workplan_data
    overlapping_pairs = find_overlapping_events(filtered_data, threshold=10)
    for pair in overlapping_pairs:
        print(f"Overlapping Pair: {pair[0][0]} and {pair[1][0]}, Start Time: {pair[0][1]}, End Time: {pair[0][2]}, Start Time: {pair[1][1]}, End Time: {pair[1][2]}")
        
    adjusted_data = adjust_event_intervals(filtered_data)
    
    categories = ['GP', 'OTHER', 'TOO-NOM', 'TOO-EX']
    category_indices = {cat: i for i, cat in enumerate(categories)}
    
    fig, axs = plt.subplots(3, 1, figsize=(10, 6), sharex=True)
    
    # Plot before filtering
    axs[0].set_title("Workplan Timeline Before Filtering")
    for item in workplan_data:
        start = item['start']
        end = item['end']
        category = item['category']
        y_pos = category_indices.get(category, -1)
        if y_pos != -1:
            # axs[0].plot([start, end], [y_pos, y_pos], marker='o', label=category if y_pos == 0 else "")
            axs[0].plot([start, end], [y_pos, y_pos], linewidth=4, label=category if y_pos == 0 else "")
    axs[0].set_yticks(range(len(categories)))
    axs[0].set_yticklabels(categories)
    axs[0].grid(True)
    
    # Plot after filtering
    axs[1].set_title("Workplan Timeline After Filtering")
    for item in filtered_data:
        start = item['start']
        end = item['end']
        category = item['category']
        y_pos = category_indices.get(category, -1)
        if y_pos != -1:
            # axs[1].plot([start, end], [y_pos, y_pos], marker='o', label=category if y_pos == 0 else "")
            axs[1].plot([start, end], [y_pos, y_pos], linewidth=4, label=category if y_pos == 0 else "")
    axs[1].set_yticks(range(len(categories)))
    axs[1].set_yticklabels(categories)
    axs[1].grid(True)
    
    axs[2].set_title("Workplan Timeline After Adjusted")
    for item in adjusted_data:
        start = item['start']
        end = item['end']
        category = item['category']
        y_pos = category_indices.get(category, -1)
        if y_pos != -1:
            # axs[1].plot([start, end], [y_pos, y_pos], marker='o', label=category if y_pos == 0 else "")
            axs[2].plot([start, end], [y_pos, y_pos], linewidth=4, label=category if y_pos == 0 else "")
    axs[2].set_yticks(range(len(categories)))
    axs[2].set_yticklabels(categories)
    axs[2].grid(True)
    
    for ax in axs:
        ax.xaxis.set_major_locator(mdates.AutoDateLocator())
        ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M'))
        # ax.legend()
    
    plt.xticks(rotation=45)
    plt.xlabel("Time")
    plt.tight_layout()
    plt.show()

def plot_workplan():
    workplan_data = getSvomWorkplan()
    filtered_data = filter_overlapping_events(workplan_data)
    adjusted_data = adjust_event_intervals(filtered_data)
    
    categories = ['GP', 'OTHER', 'TOO-NOM', 'TOO-EX']
    category_indices = {cat: i for i, cat in enumerate(categories)}
    
    fig, axs = plt.subplots(3, 1, figsize=(10, 7), sharex=True)
    
    # Plot before filtering
    axs[0].set_title("Workplan Timeline Before Filtering")
    for item in workplan_data:
        start = item['start']
        end = item['end']
        category = item['category']
        y_pos = category_indices.get(category, -1)
        if y_pos != -1:
            axs[0].plot([start, end], [y_pos, y_pos], linewidth=2, label=category if y_pos == 0 else "")
    axs[0].set_yticks(range(len(categories)))
    axs[0].set_yticklabels(categories)
    axs[0].grid(True)
    
    # Plot after filtering
    axs[1].set_title("Workplan Timeline After Filtering")
    for item in filtered_data:
        start = item['start']
        end = item['end']
        category = item['category']
        y_pos = category_indices.get(category, -1)
        if y_pos != -1:
            axs[1].plot([start, end], [y_pos, y_pos], linewidth=2, label=category if y_pos == 0 else "")
    axs[1].set_yticks(range(len(categories)))
    axs[1].set_yticklabels(categories)
    axs[1].grid(True)
    
    # Plot after adjusting intervals
    axs[2].set_title("Workplan Timeline After Adjusting Intervals")
    for item in adjusted_data:
        start = item['start']
        end = item['end']
        category = item['category']
        y_pos = category_indices.get(category, -1)
        if y_pos != -1:
            axs[2].plot([start, end], [y_pos, y_pos], linewidth=2, label=category if y_pos == 0 else "")
    axs[2].set_yticks(range(len(categories)))
    axs[2].set_yticklabels(categories)
    axs[2].grid(True)
    
    for ax in axs:
        ax.xaxis.set_major_locator(mdates.AutoDateLocator())
        ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y-%m-%d %H:%M'))
        # ax.legend()
    
    plt.xticks(rotation=45)
    plt.xlabel("Time")
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    
   plot_workplan2()